clear all

load ParEst_eh1_time0_Mex_nokids1 ParEst_eh1_time0_Mex_nokids1
load W_time1_age0_eh1_comp01.out
load W_time1_age0_eh1_comp0.out
load W_time1_eh1_comp0.out
load stock_time1_age0_eh1_comp0.out
load stock_time1_eh1_comp0.out
load trans_time1_age0_eh1_comp0.out
load trans_time1_eh1_comp0.out

Par1 = zeros(25,1);
Par1(1:24) = ParEst_eh1_time0_Mex_nokids1;
Par1(21) = 10^(Par1(21)); 
Par1(22:24) = Par1(22:24)*1000;
Par1(25) = 0;

w_f_g=W_time1_age0_eh1_comp01; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

w_mom=W_time1_age0_eh1_comp0; % Read wage moments
meanwi_1=w_mom(1); meanwf_1=w_mom(2);   meanwi_2=w_mom(3);  meanwf_2=w_mom(4);  sdwi_1=w_mom(5); sdwf_1=w_mom(6);  sdwi_2=w_mom(7);   sdwf_2=w_mom(8);
p10wi_1=w_mom(9);  p10wf_1=w_mom(10); p10wi_2=w_mom(11); p10wf_2=w_mom(12); p25wi_1=w_mom(13); p25wf_1=w_mom(14); p25wi_2=w_mom(15); p25wf_2=w_mom(16); 
p50wi_1=w_mom(17); p50wf_1=w_mom(18); p50wi_2=w_mom(19); p50wf_2=w_mom(20); p75wi_1=w_mom(21); p75wf_1=w_mom(22); p75wi_2=w_mom(23); p75wf_2=w_mom(24); 
p90wi_1=w_mom(25); p90wf_1=w_mom(26); p90wi_2=w_mom(27); p90wf_2=w_mom(28); minwi_1=w_mom(29); minwf_1=w_mom(30); minwi_2=w_mom(31); minwf_2=w_mom(32); 
maxwi_1=w_mom(33); maxwf_1=w_mom(34); maxwi_2=w_mom(35); maxwf_2=w_mom(36);
wage_moments_data=[meanwi_1; meanwf_1; meanwi_2; meanwf_2; sdwi_1; sdwf_1; sdwi_2; sdwf_2;...
    p10wi_1; p10wf_1; p10wi_2; p10wf_2; p25wi_1; p25wf_1; p25wi_2; p25wf_2;...
    p50wi_1; p50wf_1; p50wi_2; p50wf_2; p75wi_1; p75wf_1; p75wi_2; p75wf_2;...
    p90wi_1; p90wf_1; p90wi_2; p90wf_2];

Dij=trans_time1_age0_eh1_comp0'; % Read transition moments

stock_mom=stock_time1_age0_eh1_comp0; % Read stocks
stock=stock_mom';

% Transitions
Dnf_1	 =	Dij(1);
Dni_1	 =	Dij(2);
Dfn_1	 =	Dij(3);
Dff_1	 =	0;
Dfi_1	 =	Dij(4);
Din_1	 =	Dij(5);
Dii_1	 =	0;
Dif_1  =	Dij(6);
Dni_df_1 = Dij(7);    % Transitions for spouse 2 (spouse) from not working to informal when spouse 1 (head) loses job in formal sector
Dni_di_1 = Dij(8);    % Transitions for spouse 2 (spouse) from not working to informal when spouse 1 (head) loses job in informal sector

Dnf_2	 =	Dij(9);
Dni_2	 =	Dij(10);
Dfn_2	 =	Dij(11);
Dff_2	 =	0;
Dfi_2	 =	Dij(12);
Din_2	 =	Dij(13);
Dii_2	 =	0;
Dif_2  =	Dij(14);
Dni_df_2 = Dij(15);
Dni_di_2 = Dij(16);
 
HHstates_data=stock; % [FF;FI;FN;IF;NF;II;IN;NI;NN]; % joint-job state
transitions_data=[Dfn_1;Din_1;Dnf_1;Dni_1;Dfi_1;Dif_1;Dfn_2;Din_2;...
     Dnf_2;Dni_2;Dfi_2;Dif_2;Dni_df_2;Dni_di_2;Dni_df_1;Dni_di_1];

n=50;
 % solve dist of wage offers, value functions and simulate trajectory (MSM)
[Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
[Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn,iterVF]=solve_value_function_A3_C(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
[HHstates,wage_moments,transitions]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,...
    wi_1,ff_2,wf_2,fi_2,wi_2,Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn);


%Contact rates
Transitions_both = {'F-N'; 'I-N'; 'N-F'; 'N-I'; 'F-I'; 'I-F'; 'p1'; 'q1';''; 'F-N'; 'I-N'; 'N-F'; 'N-I'; 'F-I'; 'I-F'; 'p2'; 'q2'};
Estimates = [Par1(1); Par1(2); Par1(3); Par1(4); Par1(5); Par1(6); Par1(13); Par1(14);NaN;Par1(7); Par1(8); Par1(9); Par1(10); Par1(11); Par1(12); Par1(15); Par1(16)];
Transition_Rates_table = table(Transitions_both, Estimates);
writetable(Transition_Rates_table, 'Transition_Rates_table.xls')


alphaf=Par1(17);
betaf=Par1(18);
alphai=Par1(19);
betai=Par1(20);
theta=Par1(21);
a=Par1(22);
b1 = Par1(23); 
b2 = Par1(24);

HH_avg_Inc = p50wi_1; % median wage of informal men
MWP_a = a/(exp(-theta*HH_avg_Inc));
save MWP_a MWP_a

%Preference parameters
Preference = {'theta';'b1';'b2';'a';'MWP a'};
Estimates = [theta; b1; b2; a; MWP_a];
Preferences_table = table(Preference, Estimates);
writetable(Preferences_table, 'Preferences_table.xls')

% Generate mean LOG WAGES - F distribution 
mean_wf_1f=log(sum(wf_1.*ff_1));
mean_wi_1f=log(sum(wi_1.*fi_1));
mean_wf_2f=log(sum(wf_2.*ff_2));
mean_wi_2f=log(sum(wi_2.*fi_2));
mean_wage_offer=[mean_wf_1f;mean_wi_1f;NaN;mean_wf_2f;mean_wi_2f];
save mean_wage_offer mean_wage_offer

%Wage offer parameters
Wage_offer = {'alpha_formal'; 'beta_formal';'alpha_informal';'beta_informal';'';'mean_wf_1f';'mean_wi_1f';NaN;'mean_wf_2f';'mean_wi_2f'};
Estimates = [alphaf; betaf; alphai; betai;NaN;mean_wage_offer];
Wage_offer_par_table = table(Wage_offer, Estimates);
writetable(Wage_offer_par_table, 'Wage_offer_par_table.xls')


y=transitions_data;
transitions_data1=[y(3);y(4);y(1);y(5);y(2);y(6);y(13);y(14);NaN;NaN;...
    y(9);y(10);y(7);y(11);y(8);y(12);y(15);y(16)];
z=transitions;
transitions1=[z(3);z(4);z(1);z(5);z(2);z(6);z(13);z(14);NaN;NaN;...
    z(9);z(10);z(7);z(11);z(8);z(12);z(15);z(16)];

% Model fit (for paper)
Model_fit_labels = {'mff'; 'mfi'; 'mfn'; 'mif'; 'mnf'; 'mii'; 'min'; 'mni'; 'mnn';'';'';...
    'Dnf_1';'Dni_1';'Dfn_1';'Dfi_1';'Din_1';'Dif_1';'Dni_df_2';'Dni_di_2';'';'';...
    'Dnf_2';'Dni_2';'Dfn_2';'Dfi_2';'Din_2';'Dif_2';'Dni_df_1';'Dni_di_1'};
ModelFit_actual=[HHstates_data;NaN;NaN;transitions_data1];
ModelFit_model=[HHstates;NaN;NaN;transitions1];
Model_fit_table = table(Model_fit_labels,ModelFit_actual,ModelFit_model);
writetable(Model_fit_table, 'Model_fit_table_paper1.xls')

x=wage_moments_data;
wage_moments_data1=[x(10);x(14);x(18);x(22);x(26);x(2);x(6);NaN;NaN;...	
x(9);x(13);x(17);x(21);x(25);x(1);x(5);NaN;NaN;...
x(12);x(16);x(20);x(24);x(28);x(4);x(8);NaN;NaN;...	
x(11);x(15);x(19);x(23);x(27);x(3);x(7)];
w=wage_moments;
wage_moments1=[w(10);w(14);w(18);w(22);w(26);w(2);w(6);NaN;NaN;... 
w(9);w(13);w(17);w(21);w(25);w(1);w(5);NaN;NaN;...
w(12);w(16);w(20);w(24);w(28);w(4);w(8);NaN;NaN;... 
w(11);w(15);w(19);w(23);w(27);w(3);w(7)];

% Model fit (for paper)
Model_fit_labels = {'p10wf_1';'p25wf_1';'p50wf_1';'p75wf_1';'p90wf_1';'meanwf_1';'sdwf_1';'';'';...
    'p10wi_1';'p25wi_1';'p50wi_1';'p75wi_1';'p90wi_1';'meanwi_1';'sdwi_1';'';'';...
    'p10wf_2';'p25wf_2';'p50wf_2';'p75wf_2';'p90wf_2';'meanwf_2';'sdwf_2';'';'';...
    'p10wi_2';'p25wi_2';'p50wi_2';'p75wi_2';'p90wi_2';'meanwi_2';'sdwi_2'};
ModelFit_actual=log(wage_moments_data1);
ModelFit_model=log(wage_moments1);
Model_fit_table = table(Model_fit_labels,ModelFit_actual,ModelFit_model);
writetable(Model_fit_table, 'Model_fit_table_paper2.xls')

